#include <iostream>
#include <string>
#include <algorithm>
#include <cmath>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>

#include <TROOT.h>
#include <TFile.h>
#include <TTree.h>
#include <TRandom.h>
#include <TSQLServer.h>
#include <TSQLResult.h>
#include <TSQLRow.h>
#include <TClonesArray.h>

#include "GeomSvc.h"
#include "TriggerRoad.h"

using namespace std;

int main(int argc, char **argv)
{
  cout << "Exporting Run: " << argv[1] << " to ROOT file: " << argv[2] << endl;

  GeomSvc *p_geomSvc = GeomSvc::instance();
  p_geomSvc->init("geometry_R962");

  Int_t nHits_single;
  Int_t detectorID_single[50];
  Int_t elementID_single[50];
  Double_t x_single[50];
  Double_t y_single[50];
  
  TFile* dataFile1 = new TFile(argv[1], "READ");
  TTree* dataTree1 = (TTree*)dataFile1->Get("save");

  dataTree1->SetBranchAddress("nHits", &nHits_single);
  dataTree1->SetBranchAddress("detectorID", detectorID_single); 
  dataTree1->SetBranchAddress("elementID", elementID_single); 
  dataTree1->SetBranchAddress("x", x_single); 
  dataTree1->SetBranchAddress("y", y_single); 

  TFile* dataFile2 = new TFile(argv[2], "READ");
  TTree* dataTree2 = (TTree*)dataFile2->Get("save");

  dataTree2->SetBranchAddress("nHits", &nHits_single);
  dataTree2->SetBranchAddress("detectorID", detectorID_single); 
  dataTree2->SetBranchAddress("elementID", elementID_single); 
  dataTree2->SetBranchAddress("x", x_single); 
  dataTree2->SetBranchAddress("y", y_single); 

  Int_t nHits;
  Int_t detectorID[5000];
  Int_t elementID[5000];
  Double_t x[5000];
  Double_t y[5000];
 
  Int_t nMuons;
  Int_t nElectrons;

  TFile *saveFile = new TFile(argv[3], "recreate");
  TTree *saveTree = new TTree("save", "save");
 
  saveTree->Branch("nHits", &nHits, "nHits/I");
  saveTree->Branch("detectorID", detectorID, "detectorID[nHits]/I");
  saveTree->Branch("elementID", elementID, "elementID[nHits]/I");
  saveTree->Branch("x", x, "x[nHits]/D");
  saveTree->Branch("y", y, "y[nHits]/D");
  saveTree->Branch("nMuons", &nMuons, "nMuons/I");
  saveTree->Branch("nElectrons", &nElectrons, "nElectrons/I");

  Int_t nEntries = atoi(argv[4]);
  cout << "Totally " << nEntries << " events in this run" << endl;
 
  Int_t seed = 0; 
  if(argc > 5) seed = atoi(argv[5]);
  TRandom rnd(seed);
  for(Int_t i = 0; i < nEntries; i++)
    {
      if(i % 1000 == 0) cout << "Mixing event " << i << endl;
      
      nMuons = -1;
      while(nMuons < 0 || nMuons > 100) nMuons = Int_t(rnd.Gaus(44., 6.17));
      nElectrons = -1;
      while(nElectrons < 0 || nElectrons > 200) nElectrons = Int_t(rnd.Landau(24.2, 5.58));

      nHits = 0;
      for(Int_t j = 0; j < nMuons; j++)
	{
	  Int_t index = Int_t(rnd.Rndm()*dataTree1->GetEntries());
	  dataTree1->GetEntry(index);

	  for(Int_t k = 0; k < nHits_single; k++)
	    {
	      detectorID[nHits] = detectorID_single[k];
	      elementID[nHits] = elementID_single[k];
	      x[nHits] = x_single[k];
	      y[nHits] = y_single[k];
	    
	      ++nHits;
	    }
	}

      for(Int_t j = 0; j < nElectrons; j++)
	{
	  Int_t index = Int_t(rnd.Rndm()*dataTree2->GetEntries());
	  dataTree2->GetEntry(index);

	  for(Int_t k = 0; k < nHits_single; k++)
	    {
	      detectorID[nHits] = detectorID_single[k];
	      elementID[nHits] = elementID_single[k];
	      x[nHits] = x_single[k];
	      y[nHits] = y_single[k];
	    
	      ++nHits;
	    }
	}

      saveTree->Fill();
    }

  saveFile->cd();
  saveTree->Write();
  saveFile->Close();
  
  return 1;
}
